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PACS 68 . 03 . Fg - Evaporation and condensation of liquids 
PACS 44 . 35 . +c - Heat flow in multiphase systems 
PACS 61 . 46 . -w - Structure of nanoscale materials 

Abstract. - In a one-component fluid, we investigate evaporation of a small axysymmetric liquid 
droplet in the partial wetting condition on a heated wall at T ~ Q.QTc. In the dynamic van der 
Waals theory (Phys. Rev. E 75, 036304 (2007)), we take into account the latent heat transport 
from liquid to gas upon evaporation. Along the gas-liquid interface, the temperature is nearly 
equal to the equilibrium coexisting temperature away from the substrate, but it rises sharply to 
the wall temperature close to the substrate. On an isothermal substrate, evaporation takes place 
mostly on a narrow interface region near the contact line in a late stage, which is a characteristic 
feature in one-component fluids. 
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Introduction. — The wetting dynamics has been 
mostly studied for imvolatile liquids and is not well un- 
derstood for volatile liquids [1]. A well-known example 
js evaporation of a liquid droplet on a heated substrate. 
Here first-order phase transition from liquid to gas occurs 
on the interface, where latent heat is carried away from the 
interface with gas flow. For a liquid droplet in air the ra- 
dius rc{t) has been observed to decrease as rc{t) ^ {tn — t)°- 
with a = 0.5 until it vanishes at time to [2-5]. Detailed in- 
formation of evaporation and contact-line motion has been 
provided by various theoretical approaches [6-12]. In near- 
critical one-component fluids, in particular, a bubble in 
liquid was observed to be attracted to a heated wall even 
when the wall was wetted by liquid in equilibrium [13]. 
In one-component fluids, the contact angle decreases from 
zero to an apparent finite value in heat flux. 

In this Letter, we numerically study evaporation of a 
droplet in one-component fluids in the axisymmetric ge- 
ometry. As an efficient numerical method, we use the dy- 
namic van der Waals model [14], which is a phase fleld 
model of fluids with inhomogeneous temperature. The 
pressure p outside a droplet (or bubble) is nearly ho- 
mogeneous so that the interface temperature should be 
close to the equilibrium coexisting (saturation) temper- 
ature Tcx(p) even in heat flux [14,15]. Thus, near the 
contact line on a heated or cooled wall, a steep temper- 
ature variation and a large heat flux should appear, as 
theoretically studied by Nikolayev et al. [11, 12] and as 



measured by Hohmann and Stephan [16]. Therefore, the 
hydrodynamics is singular around the contact line in heat 
flux. This aspect has not yet been well studied in the liter- 
ature. On the other hand, in multi-component fluids, the 
interface temperature changes on the scale of the droplet 
size and evaporation should take place all over the surface. 

Dynamic van der Waals model. We examine 
the gas-liquid phase transition in nonstationary, inhomo- 
geneous temperature T{r,t). We start with the entropy 
functional Sf, dependent on the number density n(r, t) and 
the internal energy density e(r, f) and introduce T{r, t) 
by the functional derivative l/T — {dSb/dn)e- We assume 
that Sb = J drS is the space integral of the entropy den- 
sity S consisting of regular and gradient parts as 



S = ns{n,e)- -C\Vnf 



(1) 



where s is the entropy per particle dependent on n and e. 

The second term is the gradient entropy density, which is 
negative and is imortant near the interface. The coefficient 
C is assumed to be a constant. In the van der Waals 
theory, fluids are characterized by the molecular volume 
vq and the attractive pair interaction energy e and the 
entropy s is written as 

s = kB ln[(e/n + evonf^'^{l/vQn - 1)] + sq, (2) 

where so/ks = ln[uo(m/37r?i^)^/^] -|- 5/2 with m being the 
molecular mass. In this Letter, we assume that the fluid 
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internal energy £{, is the space integral of e. Then we have 
the usual relation 1/T — n{ds/de)n, which yields the well- 
known expression e = 3nkBT/2 — tv^TT? . More generally, 
we may assume the form = / (ir[e + iir|Vnp/2], where 
the second term represents the gradient energy density 
[14]. In this Letter, we set K = Q for simplicity. 

We set up the hydrodynamic equations from the princi- 
ple of positive entropy production in nonequilibrium. No 
gravity is assumed. The mass density p — mn obeys the 
continuity equation. 



d_ 
di 



P ■ 



-V • {pv). 



(3) 



The momentum density pv and the total energy density 
ct = e + pv^ /2 are governed by [17] 

d_ 
d 



-pv 



-V- {pw + n- a), 



(4) 



(5) 



ceiling wall 




substrate 



Here n — {^ij} is the reversible stress tensor consisting 



of the van der Waals pressure p = kBTn/(l - 
and the gradient contribution as 



Hij = p5u - CT 



{n\/'^n+\\\/n\^)5., 



(6) 



The a = {a,j} = r]{V,Vj + VjV,) + (C - 2ri/3)V ■ vS,j is 
the viscous stress tensor in terms of the shear viscosity 77 
and the bulk viscosity (. Hereafter Vi d/dxi with Xi 
representing x, y, and z. The A in eq. (5) is the thermal 
conductivity. 

We note that H^- satisfies Ylj^j^^ijl'^) — — 
eViT~^, where p is the generalized chemical potential in- 
cluding the gradient contribution. If C is assumed to be 
a constant, p, reads 



p = T{5Sb/5n)e ^ p- TCV^n, 



(7) 



where Sb — J drS is the entropy functional and p = 
p(T, n) is the usual chemical potential. In equilibrium 
the stress balance J^j ^j^ij = is equivalent to the ho- 
mogeneity of p. The interface profile n — n(x) is obtained 
from p{T, n) — CTcPn/ dx^ ~ Pc^{T), where p^^iT) is the 
chemical potential in two phase coexistence. This inter- 
face equation was derived by van der Waals [18]. 

In our simulation, however, we solved the equation for 
the entropy density S in cq. (1), 



dS_ 

'di 



(8) 



together with the continuity and momentum equations (3) 
and (4) . The right hand side of eq. (8) is the nonnegative- 
definite entropy production rate with 



a,,V,v„ eg = X{VTf/T. 



Fig. 1: Axisymmetric liquid droplet (in black) on the substrate 
in gas (in gray) for $1 = —0.1 and T = 0.875Tc in equilibrium 
in a cylindrical cell {Q < z < H and < r < L). Here the 
region r < L/2 is shown. A thermalUy insulating side wall is 
at r = L. In this work H = AOOAx = 200£ and L = 700Aa; = 
350^. 

In the literature on simulations of two phase fluids [19- 
21], it is well-known that a small velocity field remains 
nonvanishing around the interface of a droplet after long 
times even without heat input from the boundaries. It 
is an artificial parasitic fiow and its magnitude depends 
on the discretization method. In our previous work [14], 
we integrated the energy equation (5) and fig. 3 there is 
affected by such a parasitic flow. Instead, if we use the 
entropy equation (8), the entropy production rate tends 
to zero or WiVj and V^T — > at long times if there 
is no heat flow from outside. Thus, with our new method, 
equilibrium can be reached around the interface region. 

In passing, sq in eq. (2) disappears in eq. (8) owing to 
the continuity equation (4), so It is an arbitrary constant. 



Numerical method. 

with < z < H and < r = x 
a fluid. The system volume ttL^H is thus flxed 



We assume a cylindrical cell 
= ^a;2 + y2 <; ^ filled with 
In this 



axisymmetric geometry, we assume that all the variables 
depend only on z, r, and t. We integrated the dynamic 
equations on a two-dimensional lattice with H — 400Aa; 
and L = 700Aa;, where Ax = £/2 is the mesh size of the 
integration with 



£^{C/2kBVoy/^ 



(10) 



in terms of C in eq. (1) and vq in the van der Waals the- 



(9) ory. For C 



5 /3 

ksVQ we have t 



1/3 



but C remains 



arbitrary in our theory. The interface width in the follow- 
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ing simulations is then of order 6 Ax with Ax = £/2 and 
we obtained almost the same results even for Ax — £/A. 
We assumed the linear density-dependence of the trans- 
port coefficients as 77 = <^ = v^^mn and A = kgiyon, though 
they are crude approximations [22] . The kinematic viscos- 
ity 1^0 — rj/p is assumed to be a constant. Space and time 
will be measured in units of £ and 



To = f/i^o = C/2kBVoVo, 



(11) 



respectively. Away from the criticality, the thermal dif- 
fusivity Dt = X/Cp is of order uq (where Cp is the iso- 
baric specific heat per unit volume) or the Prandtl num- 
ber Pr = VqI Dt is of order unity. Hence tq is the thermal 
relaxation time on the scale of £. There arises a dimen- 
sionless number a = v^m/e£'^ and we set a — 0.06 in this 
work. Then i/q = {Qme/mf/^e or C = /csWo^^m/O.OSe. 
The velocity field v vanishes on all the boundaries. We 
control the boundary temperatures at z = and H , writ- 
ten as Tq and T^, while the side wall at r L is thermally 
insulating or {dT/dr)r=L = 0. On the substrate z = we 
imposed the boundary condition [8], 



dn 



(12) 



where $1 arises from the short-range interaction between 
the fluid and the wall assumed to be smooth [1]. We set 
dn/dz — at z ~ H and dn/dr — —O.I/vqI dX r — L. 

We first placed a hemisphere liquid droplet with radius 
100 at T = 0.875Tc, where the liquid and gas densities 
were those on the coexistence curve given by = 0.580/fo 
and rig = 0.122/t)o, respectively. We then waited a time 
interval of 4200 to let the system relax to the true equi- 
librium, where T — Q.SlbTc throughout the system and 
V — Q even around the interface, and the pressure dif- 
ference between the two phases satisfied the Laplace law. 
The liquid and gas densities attained are 0.581/wo and 
0.123/z;o, respectively, which are slightly different from the 
initial values due to the surface tension effect. In fig. 1, 
we show such a state with <I>i = —0.1, where the contact 
angle is 135°. AtT — 0.875rc, the contact angle vanishes 
for <i>i = 0.12 and the wall is completely wetted by liq- 
uid for larger $1. We determine the interface position by 
n(r, z) = (n£ -|- ng)/2. As z — > we obtain the position of 
the contact line r = rc. 

Heating a liquid droplet. — We prepared an equi- 
librium state with $1 — 0.05 in the partial wetting condi- 
tion, where the contact angle is 57°. We then raised the 
substrate temperature Tg to O.OlTc, while the temperature 
Th at the ceiling was unchanged from 0.875Tc. We set 
t = at this temperature change. The <i>i was kept fixed. 
The contact angle increased up to 67° around t — 1800 in 
accord with the experiment [13], but it slowly decreased 
afterwards being equal to 60° at t = 10800. (These two 
times correspond to the beginning and middle of an late 
stage of evaporation in fig. 4 below.) The droplet assumed 
a cap-like shape until it disappeared. 
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Fig. 2: Droplet shape (top) and velocity field (middle and bot- 
tom) for $1 — 0.05 at t = 100 and 5800 on a heated sub- 
strate with the darker regions representing liquid. Evapora- 
tion strongly takes place around the contact line, as can be 
seen from the velocity field (arrows). 



In fig. 2, the droplet shape and the velocity field are 
shown at i = 100 in an early stage and at i = 5800 in 
a late stage. Evaporation is taking place strongly in the 
vicinity of the contact line r = rc(t). In the upper panel 
of fig. 3, we show the heat fiux on the substrate. 



dT 
dz 



(13) 



z=0 



at these two times. It exhibits a sharp peak at r = rc{t), 
where the wall supplies excess heat needed for evapora- 
tion. This behavior is in accord with the theoretical result 
by Nikolayev et al. [11] In the lower panel of fig. 3 at 
t — 5800, the temperature T is nearly constant along the 
interface, as in our previous simulation [14]. It is nearly 
equal to the coexistence temperature Tcx(p), where p is 
the pressure homogeneous outside the interface region, as 
anlytic calculations demonstrated [15]. 
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Fig. 3: Upper plate: Heat flux Qbi'i", t) on a heated substrate in 
units of e£/voTo at t = 100 and 5800 with a sharp max;imum at 
the contact Une r = rdt). Lower plate: Profile of T{r, z, t)/Tc 
at t = 5800. It sharply dropls to 0.8815 near the contact line 
(r = rc{t) and z = 0), where evaporation occurs strongly. 

In our simulation. T sharply drops from Tq to Tf.x{p) 
near the contact line on the scale of the interface thick- 
ness ^, while T gradually decreases outside the contact 
line region |r — rd 3> ^. We confirmed this results even for 
other To in the partial wetting condition. The heat flux 
near the contact line is then given by Qcon = ^t^T/S^, 
where AT — Tq ~ Tex and is the thermal conductivity 
of liquid. This heat flux may be equated with the convec- 
tive latent heat flux ngT^^l^sVg in the gas region, where 
Ug is the gas density, As is the entropy difference per par- 
ticle, and Vg is the gas velocity near the contact Ine. The 
convection dominates over the thermal diffusion in this 
region. Therefore, 

Vg ~ A^AT/(en3TcxAs) 

~ voneAT/^UgT^. (14) 

In the second line we set = /cb^'o^^ and As ^ ks- This 
estimation is consistent with the numerical values of Vg 
in fig. 2. We also changed AT and confirmed the linear 
relationship Vg oc AT in the partial wetting condition. 

In fig. 4, we display the radius of the contact line rc{t) 
versus t. In the early stage t <10'^, ^c{t) decreases rapidly, 
where evaporation takes place strongly all over the surface 
(see fig. 5 below). In the late stage t >10^, it decreases 




txio^ 

Fig. 4: Time evolution of the contact line radius rc{t) of an 
evaporating droplet. For t >W^ it may be fitted to eq. (15). 

For t<lO^ it decreases rapidly due to enhanced evaporation, 
algebraically as 

r^t) a {to - ^)°•^^ (15) 

until it disappears at to — 2 x 10^. The exponent 0.42 
in eq. (15) is smaller than the exponent 0.5 for macro- 
scopic droplets in air [2 4] . The droplet volume decreased 
as {to — t)^'^"^ and the average droplet density n{t) slowly 
decreased roughly as n — Ug oc {to — t)^'^^. The latter is 
because the droplet interior is gradually heated, as can 
be seen from the inhomogeneous temperature profile in 
fig. 3. The interior density approaches the gas density be- 
fore its disappearance. The droplet interior thus changes 
in a complicated manner and we cannot present clear ex- 
planations of the exponents 0.42, 1.12, and 0.08 at present. 

Evaporation rate. — We next examine the evapora- 
tion rate. To this end we introduce the mass flux through 
the interface, 

J=n{v- Vint) ■ 1^, (16) 

where f = — |Vn|~^Vn is the normal unit vector at the 
surface from liriuid to gas and Vint ■ i^{= drc{t)/dt) is the 
interface velocity. If J is regarded as a function of the 
coordinate along the normal direction i/, it is continuous 
through the interface from the mass conservation, while n 
and V ■ V change discontinuously. Thus we may determine 
J = J{r, t) on the interface as a function of r at each time. 
In the thin interface limit J is related to the discontinuity 
of the heat flux —A • VT as 

TAsJ - u ■ [(AVT)gas - (AVT)uq] = 0, (17) 

from the energy conservation at the interface, where the 
subscript gas (liq) denotes the value in the gas (liquid) side 
close to the interface. The total evaporation rate W^^^ is 
the surface integral of J. The surface area in the range 
[r,r + dr] is 2Trdrr/ sinO, where 6 is the angle between i/ 
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Fig. 5: Evaporation rate on all the interface VKtot in eq. (18) 
and that Wcon in the vicinity of the contact line rc — 6£ < r < Vc 
in eq. (19) in units of £^/voto vs t until disappearance of the 
droplet. Evaporation occurs all over the interface in an early 
stage after heating the substrate, but it mostly occurs at the 
contact lino in a late stage. Inset: evaporation rates in the 
early stage, where sound propagation in liquid gives rise to 
oscillatory relaxation. 



and the r axis. Thus, 



tot 



2tt 



pre 

/ drrJ/: 
Jo 



sin^. 



(18) 



The particle number within the droplet iV^ decreases in 
time as dNd/dt = —W^^^. The droplet volume Vd is related 
to Nd by Vd ~ Vd/n in terms of the average droplet density 
n and is proportional to for not thin droplets. We also 
define the evaporation rate near the contact line. 



= 2lT f 

•I rc — Vu 



drrJ/smO. 



Here we set r^^ = 6£, which is twice longer than the inter- 
face width and remains shorter than before the droplet 
disappearance. In fig. 5, we show W^tot(^) ^^'^ ^con(0 
versus t. The two curves nearly coincide for t >10^. This 
demonstrates that evaporation occurs only near the con- 
tact line at long times. Also in fig. 3, along the interface 
far from the substrate z 3> ^, we can see v ■ VT = and 
recognize J = from eq. (17). In terms of Vg in eq. (14) 
we then obtain 



(20) 



If the droplet density is treated as a constant (= n^), 
eqs. (14) and (21) yield 



dt 



(21) 



which is consistent with the curve in fig. 4 at long times 
t >10^. The agreement with the decay behavior (15) be- 
comes better if we account for the slow decrease of the the 
droplet density (see the discussion below eq. (15)). 
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Fig. 6: Lcidcnfrost effect for To = 0.975Tc and Th = 0.8751;. 
Top plates: a gas layer appears between the droplet and the 
substrate (left) and the droplet is detached from the substrate 
(middle and right). Bottom plate: Profile of T{r,z,t)/Tc at 
t = 600 after droplet detachment, where T is nearly constant 
along the interface. The gas layer between the droplet and the 
substrate has a large temperature gradient. 



(19) As shown in the inset of fig. 5, W^^^ is much enhanced 
exhibiting oscillatory behavior on the acoustic time scale 
(^the droplet radius divided by the liquid sound velocity) , 
where evaporation takes place all over the interface. The 
sound velocity is 5.43 in liquid and 2.67 in gas in units 
of e/TQ (with a = 0.06) at T = 0.875Tc, so the acoutic 
disturbances propagate faster in liquid. Before the first 
peak of W^tot, a compression sound pulse emitted from the 
substrate adiabatically warms the liquid [23,24], resulting 
in a rapid rise of evaporation on all the interface. Its sub- 
sequent sharp drop at t ~ 7 then occurs with propagation 
of an expansion sound wave emitted from the contact line, 
where evaporation suddenly starts and cools the surround- 
ing liquid. These acoustic waves are mixed after the first 
minimum. 



Leidenfrost effect. — In fig. 6, we show an extreme 

case of heating the substrate from 0.875Tc to 0.975Tc at 
t = 0. Here, at the pressure in the gas, the substrate tem- 
perature is much above Tcx(p) and the thermal diffusion 
is relatively slow near the substrate since To = Tc- On a 
time scale of order 100, a gas layer emerges on the sub- 
strate and it supports a large temperature gradient. Here 
the whole liquid layer adjacent to the heater turns simul- 
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taneously to the gas. The droplet is then detached from 
the substrate (the Leidenfrost effect). Without gravity 
in this case, the ch'oplct slowly moves towards the cooler 
boundary and the droplet interior becomes cooler than the 
surrounding gas region due to evaporation. It eventually 
collides with the cooler boundary to form a thickened wet- 
ting layer (see fig. 7 of our previous paper [14]). In grav- 
ity, a droplet may be suspended in gas in a steady state. 
On the other hand, a gas bubble is attacted to a warmer 
boundary [13,14]. The lower panel of fig. 6 demonstrates 
that the temperature is nearly homogeneous all along the 
interface at Tcx(p) with disappearance of the contact line. 
Here the temperature in the gas layer between the droplet 
and the substrate steeply changes from Tq to T^xip)- Thus 
the layer is strongly absorbing heat from the substrate. 

Summary and concluding remarks. For one- 
component fluids we have examined evaporation of a small 
droplet on a heated smooth substrate in the partial wet- 
ting condition in the axisymmctric geometry. In the dy- 
namic van der Waals theory [14], we have integrated the 
entropy equation (8) together with the continuity and mo- 
mentum equations to avoid parasitic flow around the in- 
terface. Our system length is of the order of several ten 
nanometers if the mesh length Ax is a few A. 

In our simulation, the temperature exhibits a sharp drop 
near the contact line, leading to a large temperature gra- 
dient and a large heat flux localized near the contact line. 
As is evident in fig. 5, evaporation in one-component fluids 
occurs only near the contact line at long times. Here we 
should note that we have assumed the isothermal bound- 
ary condition on the substrate. For finite thermal conduc- 
tivity of the wall, however, the substrate temperature is 
lowered at the contact line and the temperature drop to 
Tcx{p) in the fiuid should become much more gradual near 
the contact line [11,12,16]. 

Phenomenologically, the evaporation rate of a thin liq- 
uid droplet in air has been assumed to be of the form [2-4] , 

Jir,t) = Jo/Vrcity-r\ (22) 

Here Jq is a constant, but its expression in terms of the 
physical parameters remains unknown. The total evapo- 
ration rate W^^^., which is the surface integral of J{r,t), 
is proportional to Tc as in our case in eq. (20). If the 
temporal decrease of the liquid density is negiec;te(l, the 
above J{r,t) yields rdt) oc {to — t)^'^ in agreement with 
the experiments [2-4]. This decay law is analogous to 
that in eq. (21) for one-component fiuids, although the 
forms of J{r,t) in the two cases are very different. In 
multi-component fluids, a velocity field induced by the 
Marangoni effect should serves to realize evaporation on 
all the interfcae [9, 10]. 
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